Minimization of a univariate polynomial

Contributed by: Benoît Legat

using DynamicPolynomials
using SumOfSquares
import CSDP

Consider the problem of finding both the minimum value of p = x^4 - 4x^3 - 2x^2 + 12x + 3 as well as its minimizers.

We can use SumOfSquares.jl to find such these values as follows. We first define the polynomial using DynamicPolynomials.

@polyvar x
p = x^4 - 4x^3 - 2x^2 + 12x + 3

\[ 3 + 12x - 2x^{2} - 4x^{3} + x^{4} \]

Secondly, we create a Sum-of-Squares program searching for the maximal lower bound σ of the polynomial.

model = SOSModel(CSDP.Optimizer)
@variable(model, σ)
@constraint(model, cref, p >= σ)
@objective(model, Max, σ)

\[ σ \]

Thirdly, solve the program and find σ = -6 as lower bound:

optimize!(model)
solution_summary(model)

solution_summary(; result = 1, verbose = false)
├ solver_name          : CSDP
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ └ raw_status         : Problem solved to optimality.
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ ├ objective_value      : -6.00000e+00
│ └ dual_objective_value : -6.00000e+00
└ Work counters
  └ solve_time (sec)   : 1.21808e-03

We can look at the certificate that σ = -6 is a lower bound:

sos_dec = sos_decomposition(cref, 1e-4)

(-3.000000004964362 - 1.999999998797403*x + 0.9999999995176428*x^2)^2

Indeed, p + 6 = (x^2 - 2x - 3)^2 so p ≥ -6.

Extraction of minimizers

We can now find the minimizers from the moment matrix:

ν = moment_matrix(cref)
ν.Q

3×3 SymMatrix{Float64}:
 1.0        0.0669168   3.13383
 0.0669168  3.13383     6.46842
 3.13383    6.46842    22.3383

This matrix is the convex combination of the moment matrices corresponding to two atomic measures at -1 and 3 which allows us to conclude that -1 and 3 are global minimizers.

η = atomic_measure(ν, 1e-4)
minimizers = [η.atoms[1].center; η.atoms[2].center]

2-element Vector{Float64}:
 -1.0000001975826238
  2.9999999542477047

Below are more details on what we mean by convex combination. The moment matrix of the atomic measure at the first minimizer is:

η1 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[1])), ν.basis.monomials)
η1.Q

3×3 SymMatrix{Float64}:
  1.0  -1.0   1.0
 -1.0   1.0  -1.0
  1.0  -1.0   1.0

The moment matrix of the atomic measure at the second minimizer is:

η2 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[2])), ν.basis.monomials)
η2.Q

3×3 SymMatrix{Float64}:
 1.0   3.0   9.0
 3.0   9.0  27.0
 9.0  27.0  81.0

And the moment matrix is the convex combination of both:

Q12 = η1.Q * η.atoms[1].weight + η2.Q * η.atoms[2].weight

3×3 Matrix{Float64}:
 1.0        0.0669169   3.13383
 0.0669169  3.13383     6.46842
 3.13383    6.46842    22.3383

Another way to see this (by linearity of the expectation) is that ν is the moment matrix of the convex combination of the two atomic measures.

Changing the polynomial basis

The monomial basis used by default can leave a problem quite ill-conditioned for the solver. Let's try to use another basis instead:

model = SOSModel(CSDP.Optimizer)
@variable(model, σ)
@constraint(model, cheby_cref, p >= σ, basis = ChebyshevBasisFirstKind)
@objective(model, Max, σ)
optimize!(model)
solution_summary(model)

solution_summary(; result = 1, verbose = false)
├ solver_name          : CSDP
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ └ raw_status         : Problem solved to optimality.
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ ├ objective_value      : -6.00000e+00
│ └ dual_objective_value : -6.00000e+00
└ Work counters
  └ solve_time (sec)   : 1.46699e-03

Although the gram matrix in the monomial basis:

g = gram_matrix(cref)
@show g.basis
g.Q

3×3 SymMatrix{Float64}:
  9.0   6.0  -3.0
  6.0   4.0  -2.0
 -3.0  -2.0   1.0

looks different from the gram matrix in the Chebyshev basis:

cheby_g = gram_matrix(cheby_cref)
@show cheby_g.basis
cheby_g.Q

3×3 SymMatrix{Float64}:
  6.25   5.0  -1.25
  5.0    4.0  -1.0
 -1.25  -1.0   0.25

they both yields the same Sum-of-Squares decomposition:

cheby_sos_dec = sos_decomposition(cheby_cref, 1e-4)

(-3.000000002974674 - 1.9999999995036628*x + 0.9999999997868867*x^2)^2

The gram matrix in the Chebyshev basis can be understood as follows. To express the polynomial $-x^2 + 2x + 3$ in the Chebyshev basis, we start by substituting $x$ into $\cos(\theta)$ to obtain $-\cos(\theta)^2 + 2\cos(\theta) + 3$. We now express it as a combination of $\cos(n\theta)$ for $n = 0, 1, 2$: $-(2\cos(\theta) - 1) /2 + 2 \cos(\theta) + 5/2.$ Therefore, the coefficients in the Chebyshev basis is:

cheby_coefs = [5/2, 2, -1/2]

3-element Vector{Float64}:
  2.5
  2.0
 -0.5

We can indeed observe that we obtain the same matrix as cheby_g.Q

cheby_coefs * cheby_coefs'

3×3 Matrix{Float64}:
  6.25   5.0  -1.25
  5.0    4.0  -1.0
 -1.25  -1.0   0.25

This page was generated using Literate.jl.